function [lb,ub] = fsel_linear_big(X,Y,A,set0,MXX,MXY,kmax)
% Forward selection for inference for linear function of coefficients
% Last edit: 3/28/17 CBH
% [lb,ub] = fsel_linear(X,Y,A,set0,MXX,MXY,K)
% X - design matrix
% Y - outcome
% A - linear combination of form A'beta
% set0 - set of initial regressors
% MXX - X'*X
% MXY - X'*Y
% kmax - number of forward selection steps to take

p = size(X,2);
varind = (1:p)';

supp_lb = set0;
supp_ub = set0;

lb = zeros(kmax,1);
ub = zeros(kmax,1);

for jj = 1:kmax
    varind_lb = setdiff(varind,supp_lb);
    varind_ub = setdiff(varind,supp_ub);
    
    ci_lb = zeros(numel(varind_lb),1);
    ci_ub = zeros(numel(varind_ub),1);
    
    for kk = 1:numel(varind_lb)

        use_lb = [supp_lb;varind_lb(kk)];
        zlb = X(:,use_lb);
%         b_lb = zlb\Y;
        Mlb = pinv(MXX(use_lb,use_lb));
        b_lb = Mlb*(MXY(use_lb));
        e_lb = Y - zlb*b_lb;
%         [~,V_lb] = hetero_se(zlb,e_lb,inv(zlb'*zlb));
        [~,V_lb] = hetero_se(zlb,e_lb,Mlb);

        use_ub = [supp_ub;varind_ub(kk)];
        zub = X(:,use_ub);
        Mub = pinv(MXX(use_ub,use_ub));
%         b_ub = zub\Y;
        b_ub = Mub*(MXY(use_ub));
        e_ub = Y - zub*b_ub;
%         [~,V_ub] = hetero_se(zub,e_ub,inv(zub'*zub));
        [~,V_ub] = hetero_se(zub,e_ub,Mub);
        
        ci_lb(kk,1) = A(use_lb)'*b_lb - 1.96*sqrt(A(use_lb)'*V_lb*A(use_lb));
        ci_ub(kk,1) = A(use_ub)'*b_ub + 1.96*sqrt(A(use_ub)'*V_ub*A(use_ub));

% OLD CODE
%        use_lb = [supp_lb;varind_lb(kk)];
%        zlb = X(:,use_lb);
%        b_lb = zlb\Y;
%        e_lb = Y - zlb*b_lb;
%        [~,V_lb] = hetero_se(zlb,e_lb,inv(zlb'*zlb));
%        BLB = sparse(p,1);
%        VLB = sparse(p,p);
%        BLB(use_lb) = b_lb;
%        VLB(use_lb,use_lb) = V_lb;
%
%        use_ub = [supp_ub;varind_ub(kk)];
%        zub = X(:,use_ub);
%        b_ub = zub\Y;
%        e_ub = Y - zub*b_ub;
%        [~,V_ub] = hetero_se(zub,e_ub,inv(zub'*zub));
%        BUB = sparse(p,1);
%        VUB = sparse(p,p);
%        BUB(use_ub) = b_ub;
%        VUB(use_ub,use_ub) = V_ub;
%        
%        ci_lb(kk,1) = A'*BLB - 1.96*sqrt(A'*VLB*A);
%        ci_ub(kk,1) = A'*BUB + 1.96*sqrt(A'*VUB*A);
        

    end
    imin = find(ci_lb == min(ci_lb));
    imax = find(ci_ub == max(ci_ub));
    imin = imin(1);
    imax = imax(1);
    
    supp_lb = [supp_lb ; varind_lb(imin)]; %#ok<*AGROW>
    supp_ub = [supp_ub ; varind_ub(imax)];
        
    lb(jj,1) = min(ci_lb);
    ub(jj,1) = max(ci_ub);
end
